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Abstract 

We propose an analytic perturbative scheme for determining the eigenvalues of the Helmholtz 
equation, (V 2 + k 2 )ip = 0, in three dimensions with an arbitrary boundary where ifi satisfies either 
the Dirichlet boundary condition (ip = on the boundary) or the Neumann boundary condition 
(the normal gradient of ip, is vanishing on the boundary). Although numerous works are 
available in the literature for arbitrary boundaries in two dimensions, to best of our knowledge the 
formulation in three dimensions is proposed for the first time. In this novel prescription, we have 
expanded the arbitrary boundary in terms of spherical harmonics about an equivalent sphere and 
obtained perturbative closed-form solutions at each order for the problem in terms of corrections to 
the equivalent spherical boundary for both the boundary conditions. This formulation is in parallel 
with the standard time-independent Rayleigh-Schrodinger perturbation theory. The efficiency of 
the method is tested by comparing the perturbative values against the numerically calculated 
eigenvalues for spheroidal, superegg and superquadric shaped boundaries. It is shown that this 
perturbation works quite well even for wide departure from spherical shape and for higher excited 
states too. We believe this formulation would find applications in the field of quantum dots and 
acoustical cavities. 
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I. INTRODUCTION 



A recent experiment has reported that the geometry of nanoscale second-phase particles 
are superellipsoids in nature A good estimation of the energyspectra of these geometries 
may help us to explore the shape of these particles more accurately. For estimation of en- 
ergyspectra of these nanoscale particles one needs to solve the Helmholtz equation in three 
dimensions (3D) with arbitrary boundaries. This is the main concern of this note and we 
believe, it is the first instance when an analytic perturbative formulation has been proposed 
to solve the Helmholtz equation for arbitrary shaped boundaries in 3D. 
The Helmholtz equation is one of the ubiquitous partial differential equations in physics 
as well as in engineering. It appears in different realms of physical and engineering prob- 
lems across diverse fields - like the study of membrane vibration, eigenanalysis in acoustic 
cavities, electromagnetic wave propagation (TE/TM modes) in a waveguide and quantum 
mechanics. Most of these problems have been pursued in order to find out the energy eigen- 
values/eigenfrequencies of the linear homogeneous Helmholtz equation subjected to different 
boundary conditions and various geometries in two or three dimensions. In the classical sce- 
nario, standard wave equation reduces to the Helmholtz equation for the sinusoidal time 
dependence. Also, in the quantum manifestation, the unitary evolution of wavefunction will 
transform the Schrodinger equation into the equation of our interest. Canonical examples of 
boundary conditions are either Dirichlet boundary condition (DBC) where the eigenfunction 
vanishes on the boundary or Neumann boundary condition (NBC) where the normal gradi- 
ent of the eigenfunction becomes zero on the boundary. Of the above mentioned problems, 
the cases of membrane vibrations, propagation of the TM modes within a waveguide and 
confinement of a quantum particle in potential well fall under the DBC umbrella and the case 
of transmission of TE modes through waveguides and acoustic cavities are the prominent 
examples of NBC. To classify the solutions of Helmholtz equation according to the nature of 
boundary condition and also as per the shape of the boundary is an intricate task. Analytic 
closed form solutions to these problems are available only for a restricted class of boundary 
geometries. The problems for rectangular, circular, elliptical and triangular (equilateral and 
right angled isosceles) boundaries are examples of such limited class in two dimensions 
For Helmholtz equation in three dimensions there are 11 such 'special' separable co-ordinate 
systems where the solution can be written as a product of three factors each of which are 
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dependent on only one co-ordinate [3]. The beautiful property of these family of separated 
solutions is that the linear combination of these solutions will eventually give us all solu- 
tions of the Helmholtz equation. Even in these restricted class most of the solutions have 
multi-parameter dependence and which makes them difficult to use. But, in most physical 
problems, one comes across boundaries for which one of more coordinates are constant not 
frozen in any of those 'special' co-ordinate systems. Explicit solutions are available in the 
literature only for spherical, rectangular and cylindrical enclosures In this context, 

finding solutions of the Helmholtz equation for any arbitrary shaped boundary becomes a 
dreadful job. 

Recently the two dimensional version of this problem has been investigated extensively in 
both the direction - analytical |6-11| as well as numerical 12-17 1 . But for the three di- 
mensional case such an analytic closed form solution by means of perturbation is not yet 
addressed. There are a few attempts made in this direction via numerical means in the 



context of chaos in general 3D billiard 18 



23 1 and also some exciting experiments to in- 



vestigate wave chaos in microwave cavities 24| and in resonant optical cavities [25j. In 
this paper we address the problem of solving the Helmholtz equations in 3D for simply 
connected arbitrary shaped boundaries. In this method, an arbitrary boundary has been 
approximated as a perturbation about a sphere of 'average radius' and expanded in terms 
of the spherical harmonics with a deformation parameter sitting in front of them. The next 
task is to improve upon the 'average radius' approximation by incorporating the correction 
terms arising due to the boundary perturbation. So, we expand the wavefunction and energy 
eigenvalues in a standard perturbative series and collect terms of different orders to get the 
respective governing equations. In a similar fashion we rewrite the boundary conditions for 
different orders. We then calculate the energy and wavefunction corrections for degenerate 
and non-degenerate states for both the cases of boundary condition separately. Finally we 
have compared our analytic results against the numerical values for spheroid, superegg and 
superquadric shapes and it has been shown that the analytic values agree quite well with 
the numerical ones for small departures from the sphere. The formulation is seen to work 
very well even for large deviations from sphere. To justify the previous statement, we have 
calculated the low lying spectra for an octahedron and a cube by perturbing a sphere and 
matched them with the numerically calculated and exactly known results respectively. 
The paper is organised as follows: The section [III prescribes the general scheme in abstract 
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sense. In section IHIl we apply these prescriptions to the case of spheroids, supereggs and 
superquadrics. A short conclusion and a few comments close the final section. 



II. FORMULATION 

The homogeneous Helmholtz equation on a three dimensional flat simply connected vol- 
ume T is given by, 

(V 2 + k 2 )ij= (V 2 + £)V = 0, (1) 

where V 2 is the Laplacian operator in 3D. We are interested in finding out the eigenspectrum 
in the interior of the bounded region for the case of Dirichlet and Neumann condition on 
the surface &T . The parameter k 2 is obtained by applying the Dirichlet boundary condition 
and can always be matched with E, the energy of a quantum particle confined in the above 
mentioned region. In this case, the function ip would be identified as the wavefunction of 
the particle. Whereas for example, in the case of acoustical enclosure one would apply the 
Neumann boundary condition on ip, which is nothing but a velocity potential and values of 
k (= y/E) would determine the pitch of the eigenmodes or the resonating frequencies. For 
the two dimensional case of a vibrating membrane it was shown by Rayleigh [2] that : 11 the 
gravest tone of a membrane, whose boundary is approximately circular, is nearly the same 
as that of a mechanically similar membrane in the form of a circle of the same mean radius 
or area." In the same spirit of the two dimensional case, we have developed a formalism to 
calculate the eigenspectrum of an arbitrary shaped enclosure by perturbing about a sphere 
of 'mean radius'. In Appendix I A 1 1 and I A 2 1 we have extended the Rayleigh's result for three 
dimensions and justify the perturbation around the 'mean radius'. 

It is natural to work in the spherical polar coordinate system (r,9,<p), where any closed 
surface satisfies the periodic condition r(8,<f) + 2tt) = r(6, <p). Now, we choose an arbitrary 
closed surface of the form r = r(6, <fi) as the boundary of our problem. The first step is to 
construct a sphere of 'mean radius' R which is given by 

IT 2lT 

R = -LJ J r(6,<i)) sin6 d9d<p. (2) 
o o 

The next step is to expand the given arbitrary shape as a perturbation over the spherical 
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surface of 'mean radius' Rq in terms of the spherical harmonics of the following form 

/ OO +<2 \ 

r{6, cf>) = Ro 1 + A E <W' <P))=Ro(l + \f(0, <P)) , (3) 



a=l b=—a 



where 



Y^eA) = \/^^]j^W pba{cose) Ya ' b{eA) = { ~ 1)b Fa&(M) ' (4) 

oo +a 
a=l b=—a 

Y£ is the spherical harmonics of order a and degree b. is the associated Legendre poly- 
nomial of order a and degree h and C^s are the expansion coefficients. Here, the spherical 
harmonics are chosen in such a manner that they are normalised to unity, i.e. 

7r 2n 




Y*(d, <f))Y c d (6, <P) sin0d0d0 = Mm- ( 6 ) 



o o 



The parameter A is a measure of the deformation, which in principle, should be small 
compared to unity in order to ensure the validation of the perturbation theory. However 
as we will observe in the following section that the perturbation works reasonably well for 
large deformation also. In the limit A — > 0, r(8,<p) = Rq, represents a sphere. We also note 
that the parameter A is an artificial one, which can be easily absorbed in the coefficient C h a . 
It is used here to keep track of the orders of perturbation. The missing term CqY® in the 
above expansion can always be absorbed in the definition of Ro given by ([3]) in agreement 
with (E]). 

Now, as per first approximation, at the zeroth order level, the eigenvalue Eq and the 
wavefunction ip Q of the system bounded by r(8,<f)) will be given by that corresponding to 
the sphere of 'mean radius' Rq 

Mr, OA) = N ntl ji{kr)P?{cos9)J m * = N U;l MY^, (j>) , (7) 

\pljRl for DBC 
E =\ n ' 1 (8) 

[alJRl for NBC 

where p = kr, N n j is a suitable normalisation constant (which is independent of m) with 
/ G N, n e N >0 , m = {—I, —(/ — !),•■■ , 0, ■ • • ,1 — 1,1} G Z and ji(p) is Z th order spherical 



Bessel function. / (= k n ^Ro) is the n th node of the I th order spherical Bessel function, 
i.e. ji((3 n ,i) = and a n j (= k n> iRo) is the n th zero of the derivative of the I th order spherical 
Bessel function, i.e. jftainj) — 0- The values of p n j and a n> i are distinct for different sets of 
n and I. There is no repetition of n values for a given value of / but m can have (21 + 1) 
values, so for I = the states are non-degenerate while for non-zero I values the states are 
(21 + l)-fold degenerate. 

Our next task is to improve upon the 'mean radius' approximation by incorporating the 
correction terms arising due to the boundary perturbation. So, we will expand if) and E as 

V> = + A^i + A 2 ^2 + , (9) 
E = E + XE l + X 2 E 2 + , (10) 

where the zeroth order terms represent the wavefunction and the eigenvalue of the unper- 
turbed state and rest of the terms are corrections of different orders to the wavefunction and 
the eigenvalue respectively. Substituting these expansions in (0Q) and collecting the coeffi- 
cients of different powers of A from both the sides we obtained the following set of equations 
after some rearrangements 



0(A°) : (V 2 + E ) iJ 

^(A 1 ): (V 2 + E )^i 

0(\ 2 ) : (V 2 + E ) iJ 2 



0, (11) 
-E^ , (12) 
-E^ x - E 2 ^ . (13) 



Eq. ffTTl) can easily be identified as the equation for a particle confined in a spherical 
box with ipo as the wavefunction corresponding to the eigenvalue Eq. Now, we will set up 
the formalism in a generic way for the degenerate and non-degenerate states for both the 
Dirichlet and the Neumann boundary conditions separately. 



A. Boundary conditions 

Like the governing equations the boundary conditions are written for different orders of 
perturbation for both the cases in the following. 
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1. Dirichlet boundary condition 
The boundary condition is 

^(r(M)) = 0. (14) 

Taylor expanding the above expression about r = i? using (j3J) and (J5J) one obtains the 
boundary condition at each order in the following form, 

O(A ): MRoA<f>) = 0, (15) 
^(A 1 ) : ^(i2o, 0, 0) + i2o/(0, 0)^ o (i? o , 0, 0) = 0, (16) 

G(A 2 ) : ^(#o, 0, 0) + Rof{0, Wi(Ro, 0, 0) + ^-f 2 (6, 0)<(/?o, 0, 0) = 0, (17) 
where 

00 +a 

/(M = ££cK(M). (is) 

a=l 6=— a 

From now on the dependence of Y£ (and /) and on the arguments (8, 0) and (r, 0, 0) 
respectively is not shown explicitly for brevity. In the above expressions, prime (/) denotes 
partial differentiation with respect to r. 



2. Neumann boundary condition 

The Neumann boundary condition is given by 

W>-n = 0, (19) 

where n is the normal at the boundary surface given by r = r(9,<f>). Expanding ( TT9|) in 
the Taylor series about r = R using (j3J) and (J5J), the conditions at different orders of A are 
obtained as follows, 

£>(A°) : Vo(^o) = 0, (20) 

^(A 1 ) : ViM + Rom(Ro) ~ ^-tto(Ro) ~ p 1 2 J te) = 0, (21) 

iio -Ho sm 

G(A 2 ) : ^2(^0) + RoM(Ro) + IrUV(Ro) - ^Ni(Ro) 

+ ^rffMRo) - — SW&M + -^^ffMRo) = o. (22) 

iio -Ro sm 6* it sm 6/ 
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In the above expressions prime ('), hat (~) and check (") denote partial differentiation with 
respect to r, 6 and <fi respectively. In the following we will discuss separately the non- 
degenerate and degenerate cases. Since, the governing equations are same for both the cases 
of DBC and NBC, the general solutions are same with differences only in the argument 
p, respective coefficients and normalisation constants which are dictated by the respective 
boundary condition. 



B. Non-degenerate states (I = 0) 

1. Dirichlet boundary condition 
For I = state we have, 

V>o = N nfi j (p)Y ° = N nfi j (p) , (23) 

where p = (/3„,o/-Ro) r , jo is the zeroth-order spherical Bessel function and N n>0 = iV ni0 F = 
{2nR^il((5 n fl))~^ is the normalisation constant. E is obtained from (jSJ) with I = and an 
appropriate value of n, as ipo is constrained by the boundary condition (115p . 
The first order correction to the wavefunction is obtained by solving ffT2"j) . Substituting the 
expression for ipo from (1231) into ffT2l . the most general solution is given by 

^ = £ E ^pW + A °o Y oMp) - ^N^Up). (24) 

p=l q=~p " 

The last term in the above solution is the particular integral to the differential equation (!T2|) 



Hi this specific 



due to the inhomogeneity. Taking a hint from the two dimensional case 

ay 

closed form of the particular solution was inspected. Imposing the boundary condition ( fi~6l) 
over ipi and matching the coefficients of different orders of spherical harmonics, we obtain the 
first order eigenvalue correction (comparing the coefficient of Y®) as well as the coefficients 
(comparing the coefficient of Y£) of the following form 

A% = p N n ,oCfj^; Po = p\ r=Ro = fa, (25) 

Ei = 0. (26) 

The remaining constant which is not required for our present purpose, can be fixed by 
normalising the corrected eigenfunction over the entire volume T. From (126]) it is evident 



that any possible non-zero correction to the eigenvalue eigenvalue will come from higher 
orders. 

The second order calculation will replicate that of the first order. Plugging the value of 
Ei = in f[T3"j) we obtained the second order correction to the wavefunction as 

oo +p 



4* = E E B lh(p)Y p q + B° Y °Jo(p) - N n , Up). (27) 

p=l q=-p " 



Implementing the boundary condition ( 1T7|) over ip2 we have estimated the second order 
correction to the eigenvalue by collecting the constant terms in flT7|) . as 



r^V^'^wi' (28) 



uq _ PoJi(po) A o v o nq NnfiPoMpo)^ ST / (2o+l)(2fe + l) nhnq -u 



a—p 
19-61 



x ( a ife00|p0)(oife6(g - 6)|pg) ( 1 + PoJ ° (p °' 

V Ja(Po) 



(29) 



where 



a—p 
|g-b| 



; Maximum (\a — p\, \q — b\) } p$ is the same as defined earlier in (1251) and 
(nin2mim2\nm) gives the Clebsch-Gordan coefficient for the decomposition of \ni,n2,n,m) 
in terms of \ni, mi) \ri2, m2).The remaining constant can be extracted out by normalising 
the corrected wavefunction up to the second order. At each order, the corrections to the 
wavefunction and eigenvalue are expressed in a generic manner in terms of the expansion 
coefficients C h a of the boundary perturbation. This formalism works for any type of closed 
surface which is slightly perturbed from a spherical one. So, for a given surface, we first 
need to calculate the 'mean radius' Rq and subsequently C^s to estimate the second order 
correction to the eigenvalue. 

2. Neumann boundary condition 
For I = state we have, 

V>o = 4,o Jo(p)Y ° = N nfi 3o (p) , (30) 

where p = (a n fl/Ro)r and N n $ = iV n) oF = {2-KRQjl(a n $)) ~2 is the normalisation constant. 
E is obtained from flB]) with I = and a given value of n. 



Solution of (fl2l) yields the first order correction to the wavefunction. Plugging ipo into (TT2I) . 
we obtained the general solution, like the DBC case 



oo +p 



fi = EE Aq Mp)y P q + A°r Q io(p) - ^-N n , ji(p). (31) 

p=l <J=-p " " 

Imposing the boundary condition (T2~T1) on ^ we have obtained 

Al = p Q N nfl Cl J j^ y p = p\ r=Ro = a nfl , (32) 
Ex = 0. (33) 

This establishes the 3D generalisation of 2D Rayleigh's theorem mentioned earlier, viz. 
"tones of symmetric modes of an arbitrary shaped cavity which is not very elongated is very 
nearly equal to that of a sphere of the same mean radius" . The remaining constant A$ can 
be fixed by normalising the corrected eigenfunction over the entire volume T ■ 
Substituting the value of E\ — in ( TTBl we obtained the second order correction to the 
wavefunction as 



oo +p 



V* = E E Hh(P)Y p q + B° Y °UP) - ^Wp). (34) 

p=l q=-p 

As ^2 is restricted by the boundary condition (1221 . we get the second order correction to 
the eigenvalue as 



^-^i- ar l 1 + l^rJ' (35) 

™ _ PoJo(Po) tivo™ j_ N n ,oPojo(po) / (2a + l)(2fc + 1) „ 6 „ 9 _ 6 

= W^ c * + ~«(ao- £ V c ' Ct 

k \\q-b\\ 

l i.nni n\/ / 1,/ mi \ Ji , Poja(po) . k(k + 1) - a(a + 1) - p(p + 1) j a (p ) \ , . 

x (afcOO p0)(afc&(g - b)\pq) < 1 + + — — — > . (36) 

I Ja(Po) 2p JaiPo)) 

By collecting the coefficients of Y£ in (122"]) we have extracted out B^. However, these con- 
stants are required only for the calculation of higher order correction terms. The remaining 
constant Bq in ip2 can be obtained from the normalisation condition of the corrected wave- 
function up to second order. We observe that these results are beautiful generalisations of 



our earlier work in two dimensions 



111 
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C. Degenerate case (I ^ 0) 



1. Dirichlet boundary condition 

For I 7^ the wavefunction of the unperturbed state is given by 

A = N ntl ji(p)Yr, (37) 

where N n j = \/2(Rljf +1 ((3 n: i))~^ is the normalisation constant and p = ((3 n) i/Ro)r. Unper- 
turbed eigenvalue Eq is dictated by (jHJ) for a given choice of I and n. All the levels with a 
given non-zero / value are (21 + l)-fold degenerate. For mathematical simplicity, we assume 
that C h a = for b ^ 0, 

oo 
a=l 

This corresponds to the boundaries with azimuthal symmetry (0 independent). So, for the 
degenerate case we will consider only the axisymmteric surface. As a result of that we can 
treat these degenerate cases in the same way as the non-degenerate case. 
The first order correction to the wavefunction can be obtained, by solving the equation (fT2l) . 

as 



^i = EE Aq MpK q + [ATM - ffr^wCp)} Yi m - (38) 



p=0 q=—p 



Now, following similar procedure as that of the non degenerate state we have using boundary 
condition ( IT5|) . 



a; = ft^to) g J W±m±l> Ck (H00| P 0> <H0m|pm>, (p + I) (39) 



fc=|J-p 



|i = _y / (4fc + l) ^ ( (2 A;)/00|/0)((2A;)/0m|/m), (40) 

where pz = /3 n> j and all other A-^'s for q ^ m are zero. The remaining constant A™ can 
be calculated from the normalisation condition. Unlike the non-degenerate states here we 
obtain a non-zero correction to the eigenvalue even at first order. 

11 



The second order correction to the wavefunction yields, 



p=0 q 



2 



Imposing boundary condition (ITT)) and collecting the constant terms we obtain second order 
eigenvalue correction as 

|i = V v/(2a+ o 1)(2&+1 -aa(a&00|fc0) 2 (A:/00|/0)(H0m|/m) 

u u a,b=lk=\a-b\ 



+ V V + l)(2fc + 1) Pdp(pi) c^Cfe <t^Z00 |p0> <7iZ0T^|prra> (fcpOQ | Z0> <fcp0ryt|Z7yi) , 

^— ' Z7T Jp(Pl) 

p=\m\ n,k=\l—p\ 

where p/ = The coefficients B^s are quite complicated and are necessary for evaluating 
the second order wavefunction and the third order eigenvalue corrections. Now usage of the 
boundary condition ( ITTj) extract out the coefficients B£. Looking at the first order degenerate 
case here we can predict that B!* is non-zero only for q = m. The constant Bp which is not 
required for the present purpose, can be obtained from the normalisation condition of the 
corrected wavefunction up to respective orders over the entire enclosed volume 7~. 



2. Neumann boundary condition 

For I 7^ the wavefunction of the unperturbed state is given by 

^ = N n ,i 3i(p)Yr, (41) 

where p = (a nj i/R )r and N n> i = {\/2a n) i/ ii{a n ^)){B^ ) {a^ l h — Z(7 + l))} - 2 is the normalisation 
constant. Eq is given by flH)) for an appropriate value of I and n. Due to the simplification 
mentioned earlier in the case of degenerate states for DBC, the boundary condition for NBC 
will reduce further as / does not contain any dependence for axisymmetric surface. So, 
we will set the terms containing / to be zero in the expression for NBC of different orders 
in (ED and (121. 



The first order correction to the wavefunction was obtained as 

oo p 



V* = E E + [aTjM - §^wG°)} Y i m - ( 42 ) 



p=0 q=—p 

p+l 
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Now, mimicking the earlier steps we have obtained using the condition (I2TT) . 



x r ? + ^ + i)-Ki + i)-P(P + i)\ (43) 

f = " E \/^^ ^ <(2t)i00|/0><(2*)/0m|im) {l + ^|^} . («) 



where = a nj ;. Like the earlier case of DBC, here also all others A^s for q ^ m are zero. 
Here also we expect a non-zero correction to the eigenvalue even at first order. 
The second order correction using ( fl3l) yields 



v» =E E (^ Jp(p) " 2i|^ Jp+i(p) ) 



p=0 q=—p 



Applying the boundary condition (122]). we have extracted out the second order eigenvalue 
correction as 

£ 2 _i 3i(i+in my 

£ 4 \p2 -/(/ + !) J^J 

/(/ + !) 1 /SA ^ /(4fc + l 



# " 1(1 + 1) 

oo 2/ 



1)} (|^) EV^^^ fc ((2A;)/00|/0)((2fc)/0m|/m) 



+ f V \/( 2q + 1)(2& + 5c a C , 6 (o600|A;0) 2 <fcf00lf0)<A;j0m|jm) 

< J < J sir 



k(k + 1) - 21(1 + 1) 
. + 2(# -/(/ + !)) 



2tt 

a, b k= 
=1 |o-6| 

oo i+p /(2n + 1V2A; 4- 1) 

V VC„C* t - - (kpQO\lO) (kpOm \lm)(nlOO\ P 0) (nlOm \pm)x 

p=|m| n,k= 

p¥=i \ 1 ~p\ 

i | rz(n + l) + Z(Z + l)-p(p+l) ) ( x | 2fl» + fc(fc + l)-p(p + l)- 1(1 + 1) Jp (p ; ) 

Few low-lying energy levels of a sphere of unit radius, subjected to either DBC or NBC, are 
tabulated in Table [fl 
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2(pf -1(1 + 1)) J I 4 ft* (p. 



Table I. Energy spectra of a unit sphere for DBC and NBC. 



DBC NBC 







Energy levels 




Eq 


States 






Energy levels 






n 


I 


771 


deg 


(= /3/ 2 J 

V " t ,71/ 




n 


I 


m 


deg 


(= «/ J 


1 








1 


9.8696 


Ground 


1 


1 


(0,±1) 


3 


4.3330 


1 


1 


(0,±1) 


3 


20.1907 


First 


1 


2 


(0,±1,±2) 


5 


11.1696 


1 


2 


(0,±1,±2) 


5 


33.2175 


Second 


1 








1 


20.1907 


2 








1 


39.4784 


Third 


1 


3 


(0,±1,±2,±3) 


7 


20.3771 


1 


3 


(0,±1,±2,±3) 


7 


48.8312 


Fourth 


1 


4 


(0,±1,±2,±3,±4) 


9 


31.8853 


2 


1 


(0,±1) 


3 


59.6795 


Fifth 


2 


1 


(0,±1) 


3 


35.2880 



III. APPLICATION TO SIMPLE CASES 



The general prescription having been sketched above, we now calculate the energy levels 
when the boundary surface is spheroidal, superegg or in general superquadric in shape. These 
are the most natural deformation of a sphere suggested by nuclear physics or experiments 
on nano scale particles Numerical results have been obtained using finite element and 
finite difference method and plots were generated using Mathematica. A direct comparison 
between the analytical results and the numerical results has been made and it shows good 
agreement for a wide range of deformation for both the boundary conditions. 



A. Spheroidal boundary 

A spheroid is a special case of an ellipsoid where two axes out of three are of the same 
length. The equation of a spheroid in Cartesian co-ordinates with z-axis as the symmetry 
axis is given by 

^ + 5 = i. («) 

a c 

where r a and r c are called equatorial radius and polar radius respectively. Now, for r c < r a 
the spheroid is called oblate, for r c > r a it is called prolate and naturally for r c = r a it 
reduces to a sphere. In spherical polar coordinates the equation of a spheroid is 
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Figure 1. Spheroids for different ratios of polar radius (r c ) to equatorial radius (r a ). Here, the 
shapes have been plotted for different values r c /r a = 0.75 (left), 1.0 (centre) and 1.25 (right). The 
figure in the centre is a sphere while the left one is an oblate and the right one is a prolate. 




r(9, 0) = T " =, (46) 

with r a > 0, 6 G [0, 7r] and G [0, 27r]. In can be easily concluded from (I46D that an oblate 
(a prolate) is a surface of revolution obtained by rotating an ellipse about its minor (major) 
axis. The shapes of spheroids for different values of r c /r a are depicted in Fig. (CQ). The 
'mean radius' Ro for a spheroid is given by 

Ro = ™< shm- 1 (VlErf) 

Since it is an axisymmetric object, it has an azimuthal symmetry. So, the expansion coef- 
ficients C„s are non-zero only for 6 = and even values of a. The magnitude of them fall 
rapidly after a few terms in the expansion. Using these expansion coefficients we have cal- 
culated the first few energy levels for spheroidal boundaries in the range 0.75 < r c /r a < 1.25 
and the comparison between the analytic perturbative values and their respective numerical 
values for DBC and NBC is shown in the Figs. (|2J) and ([3]) respectively. In these plots, an- 
alytical values are denoted by lines of different styles whereas their numerical counterparts 
are displayed by dots of different shapes. Throughout the deformation the volume of the 
spheroid remain same. 
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(47) 



Figure 2. Comparison of the eigenvalues obtained numerically (denoted by points) and analytically 
(denoted by lines) for a spheroidal boundary satisfying DBC for the first 29 states. 
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Figure 3. Comparison of the eigenvalues obtained numerically (denoted by points) and analytically 
(denoted by lines) for a spheroidal boundary satisfying NBC for the first 39 states. 
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B. Superegg boundary 



A superel 



ipse 26], a compromise between a rectangle and an ellipse, is also known as 



Lame curve 27] and is represented by the equation 

M + = (48) 
a* o* 

with t, a and b are positive numbers. A superellipsoid is a generalisation of superellipse 
into three dimensions where both horizontal sections and vertical cross-section through the 
center are superellipses with different exponent. The equation for a superellipsoid is given 
by 

(K + m n/, +\± =1 , ,49) 

with t, n, a, b and c are positive numbers. It is clear from (149]) that for a fixed value of z, 
the horizontal cross-section is a superellipse with exponent t and the vertical cross-section 
through center is also a superellipse with a different exponent n. 

A superegg [28fl is a special case of superellipsoid where a = b and t = 2. Thus, the equation 
of a superegg becomes 

aA 2 + y 2 



„ 

where a is the horizontal radius at the equator and 2c is the height of the superegg with 
exponent n. It is to be noted that, a superegg with the exponent n = 2 becomes a spheroid 
and along with that for the special case a = c, it reduces to a sphere of radius a. 
The equation of a superegg in spherical polar coordinates is given by 

^' ® = „ „ ^ ■ „, \l7n> ( 51 ) 

^ | cos I n | | smfl |nW» 

with c, a and n > 0,9 E [0, it) and G [0, 2ir]. We have considered here the case a = c = 1 
only. It is clear to understand that a superegg is a surface of revolution by rotating a 
supercircle about either of its in-plane axes. The shapes of superegg for different values 
of the exponent n are displayed in the Fig. (J3J). Since it is also an axisymmetric object, 
the expansion coefficients, C b a s are non-zero only for 6 = 0. Moreover they survive for even 
values of a and converge quickly. For the superegg shaped boundary we have calculated 
the first few energy levels using these expansion coefficients in the range 1.5 < n < 3.0 of 
the superegg exponent. The comparison between the analytic perturbative values and their 
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numerical counterparts have been illustrated in Figs. (JSJ) and fl6]) respectively for DBC and 
NBC. Here also, the analytical values are denoted by lines and the numerical ones are by 
dots. 

Figure 4. Shape of a superegg with a = c = 1 for different values of the exponent n. Here, the 
figure in the centre is a sphere of unit radius (n = 2) while the left and the right ones are supereggs 
for n = 1.5 and n = 3 respectively. 



C. Superquadrics boundary 

The equation of a general superquadric is given by 

\x\ a + \y\ b + \z\ c = 1 (52) 

where a, b and c are real positive numbers. We restrict ourselves to the case a = b = c = t 
only. For the exponent t = 1, it becomes an octahedron, for t = 2, it is a sphere and as 
t — > oo it converts to a cube. 

The equation of a superquadric in spherical polar coordinates is given by 

r(0, 0) = ; ri7? , (53) 

(| sin 6 cos 0|* + | sin 6 1 sin 0|* + | cos^l') 

with t > 0, 9 G [0,7r] and G [0, 2tt]. Our area of interest lies in the range 1 < t < 20. 
It is clear that the azimuthal symmetry is broken here unlike the above two cases. So, 
for this shape the expansion coefficients will now be non-zero for even values of a > 4 
and b mod 4 = 0. Since the general formalism is done for non-degenerate cases only and 
degenerate case requires axisymmetric perturbation, we restrict ourselves only to the non- 
degenerate cases. The shapes of the superquadric for different values of t are shown in Fig. 
([7]). For this boundary, we have calculated the eigenvalues for first three non-degenerate 
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Figure 5. Comparison of the eigenvalues obtained numerically (denoted by points) and analytically 
(denoted by lines) for a superegg boundary satisfying DBC for the first 46 states. 




1.5 2.0 2.5 3.0 

Superegg exponent (n) 

states (i.e. ground state, 9 th excited state and 45 th excited state) for DBC incorporating the 
second order eigenvalue corrections in the range 1 < t < 20 and compared them against the 
numerical results in Fig. ([8]). 
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Figure 6. Comparison of the eigenvalues obtained numerically (denoted by points) and analytically 
(denoted by lines) for a superegg boundary satisfying NBC for the first 39 states. 




1.5 2.0 2.5 3.0 



Superegg exponent (n) 

IV. CONCLUSIONS 

In this paper we have described an analytic formalism to determine the eigenspectrum of 
the Helmholtz equation in three dimensions for an arbitrary boundary subjected to either 
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Figure 7. Shape of a superquadric for different values of the exponent t. Here, the figure in the 
centre is a sphere of unit radius (t = 2) while the left one is an octahedron (t = 1) and the right 
one is a superquadric for t = 20. 




DBC or NBC. We have tested our perturbative results against the numerical values for two 
axisymmetric shapes — spheroid, superegg and superquadric geometries. This is quite an 
achievement as numerical method was the only available method to solve the Helmholtz 
equation for such boundaries. In two dimensions, the Fourier decomposition of boundary 
perturbation was an obvious choice to take care of wide variety of boundary geometries. 
In three dimensional case for a general shape, to express the boundary asymmetries on the 
surface of a sphere, the natural basis is spherical harmonics. This specific form of boundary 
perturbation in terms of spherical harmonics will be effective to address a broad class of 
boundaries. In fact, the knowledge of spherical harmonic expansion coefficients and the 
'mean radius' for a given boundary is sufficient to get the eigenvalue and the wavefunction 
corrections. Also, to write down the corrections for eigenvalue and wavefunction exactly 
in closed form at each order of perturbation is indeed extraordinary. It is evident that the 
dominant contributions are coming from the first few orders. In the first order correction 
to the eigenvalue, one C\ appears in the expressions, whereas in the second order correction 
we have a couple of C\% appearing in the expression of eigenvalue. So, the convergence of 
spherical harmonic expansion coefficients ensures the convergence of the perturbative series 
for the eigenvalue correction. The only error in the series for eigenvalue correction is due 
to the truncation of the perturbative series up to second order. The perturbative scheme 
can handle both the degenerate and non-degenerate states in a same manner with ease for 
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Figure 8. Comparison of the eigenvalues obtained numerically (denoted by points) and analytically 
(denoted by lines) for a superquadric boundary satisfying DBC for the first three non-degenerate 
states (1 = 0). 
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axisymmetric cases. Our method, even though perturbative in nature, works quite well 
for higher excited states and also for large deformation (which is against the very essence 
of the validity of a perturbative formalism). With the inclusion of higher order correction 
terms, which in principle can be obtained modulo their lengthy and cumbersome algebra, 
matching will definitely improve for the higher excited states. The (21 + 1) fold degeneracy 
of a degenerate state (with a given I value) for the spherical boundary is lifted due to the 
deformation of a sphere. For a degenerate state, the eigenvalue corrections are equal for the 
same \m\ values out of all possible m values. As a result of it, I number of levels still remain 
degenerate and we get (/ + 1) distinct energy levels. As an example we can say for / = 2, 
we should have 5 distinct energy levels but we get only three - one for / = 2, m = and 
the other two levels are for I — 2,m — ±1 and I = 2,m = ±2 respectively. For DBC, in 
the case of spheroid and superegg boundaries we have compared numerical and analytical 
energy levels up to 29 th and 46 th state (including degeneracies) respectively. Similarly for 
NBC first 39 states (including degeneracies) are compared. The matching of eigenvalues 
between analytic and numerical values are outstanding for superegg compared to spheroid 
for both the cases which is quite like their 2D analogues where the matching for supercircle 
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was better than that for ellipses js, 10|. For a larger deformation and for higher excited states 
in the case of spheroidal boundary the matching becomes less satisfactory. In these cases, 
the numerical methods may not be very accurate. It is evident from the plots that while 
the analytical method suggests exact degeneracy of some of the higher states throughout 
the entire range of deformation, the numerical counterparts split after certain range. So, 
the apparent disagreement is due to the limitation of numerical schemes rather than the 
failure of perturbative method. In case of superquadric geometries, the matching between 
the perturbative value and numerical value for the ground state is extraordinary. Maximum 
error occurs for octahedron which is ~ 5% and this is acceptable as octahedron can hardly 
be treated as a perturbation from sphere. On the other end, as t -> oo, the superquadric 
converts to a cube. So, for high value of t, say 20, the shape looks like a cube and for that 
shape we have compared analytic values against numerical ones and the discrepancy is ~ 2% 
and the eigenvalues tend towards the respective values for the cube which are |7r 2 , 47T 2 and 

In conclusion, we observe that the spherical harmonic expansion of the boundary per- 
turbation will take care of wide variety of boundaries in three dimensions. To be able to 
write down the closed form solution at each order of perturbation is a unique feature of this 
method. 
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Appendix A: Generalisation of Rayleigh's theorem in 3D 

In our formalism, the energyspectrum of an arbitrary boundary has been approximated as 
a perturbation about a sphere of 'average radius' and perturbative closed-form corrections to 
the equivalent spherical boundary at each order of perturbation have been obtained. Here, 
we will justify the above assumption for both DBC and NBC extending Rayleigh's proof for 
2D Q . 
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1. Particle in a 3D box and DBC 



In the following, we will investigate how the wavefunction of a quantum particle confined 
in a spherical box is changed by a slight deviation from the exact spherical shape. 
Independent of the shape of the boundary, ip, the wavefunction of the particle satisfies the 
wave equation (in spherical polar co-ordinate system) everywhere inside the box 

gV + 2fo/^ 1 (d 2 j) | cot0 &P | 1 d 2 VA | /c 2^,_ , Al x 

Q r 2 r Q r r 2 y QQ2 QQ g j n 2 Q Q^p, J 

where k 2 is proportional to the energy of the particle to be determined subjected to the 
boundary condition ip = on the boundary. Now, ip can always be expanded in the following 
series 

oo I 
1=1 m=—l 

where ipo, ' " , i>i are functions of r only. Plugging (1A2I) in (IA1H we obtain that, 

1=0 m=-l K v 7 ' 

which gives us ipi oc ji(kr) as the other solution of fl A3[) blows up at the origin. So, the 
general solution is given by 

oo I 

V = a jo(kr) + J2Y, a ? ^r) Y™. (A4) 

1=1 m=—l 

The boundary condition dictates us that i/j has to vanish on the boundary for all 9 and <p. 
Now, in the case of a nearly spherical well the radius vector is almost constant and can be 
approximated as r = Rq + 8r, where Sr is a small quantity which can be function of both 9 
and 0. So, the boundary condition (ip = 0) now becomes 

oo / 

a Jo (k(R + 5r)) + J2J2< ^(*(flo + 5r )) ^ = °' ( A5 ) 

1=1 m=—l 

Taylor expanding flA5j) about r = Rq we get, 

oo £ 

a MkRo) + k6r j' (kR )} + ]T ^ U(^o) + fc 5r j'^kRo)} Y t m = 0. (A6) 

Z=l m=—l 

Let us first consider the nearly symmetrical normal modes and for which we can approximate 

^ = ao ja(kr). 
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All the remaining coefficients will be small compared to a$ as the boundary is deviated a 
little from the sphere. Taking the terms up to first order of smallness, we simplify flA6j) as 

oo I 



a {jo(kRo) + kSr j' (kR )} + J2J2 a ? M^) Y ™ = °- 



(A7) 



1=1 m=—l 

Now, integrating the above expression f]A7|) with respect to 9 and <fi between the proper 
limits, we obtain the following for any arbitrary ao as 

TT 2n 



Jo(kRo) + 



k j' (kR c 
4tt 



Sr sin 9 d9 1 



0. 



0=0 4>=o 



or 



Jo 



7T 2"7T 



k\R + - 



Sr sin 9 d9 d<fi 



6=0 0=0 



Jo [k (R 0+ (5r))] = 0, 



(A8) 



where 



7r 2ir 



(Sr) 



An 



Sr sin 9 d9 d<p, 



=o 0=0 



is the average value of the change of radius from the exact spherical shape. (1A8|) shows that 
the ground state energy of a particle in a slightly deformed spherical box is nearly equal to 
that of a spherical box with radius Rq + (Sr) (= mean radius of the deformed box). 



2. Standing waves in 3D enclosures and NBC 

We now investigate the effect of slight deviation from the exact spherical shape on the 
normal modes of a spherical cavity. 

Independent of the shape of the boundary, ip which is identified as a velocity potential, 
satisfies the wave equation (lAlj) everywhere inside the sphere. Here k is proportional to 
the frequency of the normal mode to be determined subjected to the boundary condition 
d n i/j = 0, i.e. the normal component of the velocity vanishing on the boundary. Like the 
earlier case, if) can be expanded in the series given by ( ]A2j) and ()A4j) gives the appropriate 
general solution. The boundary condition dictates us that d n if> has to vanish on the boundary 
for all 9 and (p. Now, in the case of a nearly spherical cavity the radius vector is almost 
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constant and can be approximated as r = Rq + Sr, where Sr is a small quantity which can 
be function of both 9 and 0. So, the boundary condition (d n tp = 0) now becomes 



oo I 



OO I 



a-o Jo 



> (k(Ro + Sr)) + EE a ™ ilM* + 5r » Y i m ~ i E E a ? + 6r » Y ™ 5 ' 



1=1 m=—l 



1=1 m=—l 



1 



oo I 



- 2 sin 2 9 



E E a™ Ji(*(*b + ^ m ^ = 0, 



(A9) 



Z=l m=—l 



where the notations are explained in the text after ff22l) . Taylor expanding (1A9|) about 
r = Rq we get, 



oo Z 



(=1 m=—l 



a { 3o (kRo) + k Sr j' \kR )} + E E < tfi( kR o) + k 6r j'i( kR °)} *T ~ ~Y 



+k SrjKkRo)} YTSr- 



{ 3l (kR ) + k Sr jKkRo)} Yr Sr = 0. 



(A10) 



r 2 sin 2 9 

Let us first consider the nearly symmetrical normal modes and for which we can approximate 

^ = a o jo(kr). 

All the remaining coefficients will be small compared to ao as the boundary is deviated a 
little from the sphere. Taking the terms up to first order of smallness, we simplify ( lAlOj) as 



oo I 

a {f (kRo) + k Sr fo'ikRo)} + E E < M kR °) Y ™ = °" ( AU ) 

(=1 m=—l 

Now, integrating the above expression (lAllf) with respect to 9 and between the proper 
limits, we obtain the following for any arbitrary ao as 



fo(kRo) + 



k fo'jkRo 

47T 



7T 2-7T 



Sr sin 9d9 d(f> = 0, 



6»=0 <j)=0 



or 



Jo 



7r 2n 



k Ro 



47T 



Sr sin0d0d</> 



=0</>=0 



f [k(R + (5r))] = 0, 



where again 



7T 2-7T 



(Sr) = — / / Sr sin 9 d9 d(f>, 
47T J J 



=o 4>=o 



is the average value of the change of radius from the exact spherical shape. The above result 
proves that the frequencies of the vibration of symmetrical normal modes are approximately 
the same as those of a spherical enclosure with the radius having the 'mean value', Rq + (Sr). 
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Appendix B: Some useful formulae for Spherical harmonics and Clebsch-Gordan 
coefficients 



We have used the following two identities ([H [TT]) for the simplification of the product 
of spherical harmonics (the notations are explained in the text earlier). To best of our 
knowledge the identities |V] and I VI I are reported for the first time. 



I. 



\h-h\ 

[mi+ma 



II. 



V m lv m 2V m 3 _ \^ ST^ / (2/l + l)(2/ 2 + l)(2/ 3 + 1) ,,„ |+ „, , + ,„ , 

Y h Y h Y h ~ 2^ \/ 16tt 2 (2A; 2 + 1) fc2 



fci=[ l ,1_la l 

| |mi+m2| 



|fel-i 3 | 
|mi+m2+m3| 



fe 2 = 

(/i/ 2 00|/ciO)(/i/ 2 mim 2 |/ci(mi + m 2 ))(fci/ 3 00|fc 2 0)(£;iZ 3 (mi + m 2 )m 3 |fc 2 (mi + m 2 + m 3 )}; 

III. 

<9V 6 1 r "i 
^ = -gf = g |>/(a-6)(a + 6+l) F Q 6+1 - V( a + 6)(a - 6 + 1) ij" 1 } ; 

IV. 

(aa00|00)(aa6(-6)|00) = (-1) 6 (2a + l)" 1 ; 

V. a. 



\/(/ + m)(l — m + 1) (a/l(77i — l)\pm) + y(/ — m)(Z + m + l)(a/(— l)(m + l)|pm) 
a(a + 1) + (l-p)(l+p+l) 



(a/0m|pm); 



b. 



\/(/ + m)(/ - m + l)(a/l(m - l)|/m) + a/(Z - m)(/ + m + l)(a/(-l)(m + l)|/m) 
= — \fa{a + l)(a/0m|/m); 



28 



VI. 

y/b(b+l)(abO(-l)\k(-l)} + y/a(a + l)(ab(-l)0\k(-l)) = y/k(k + Tj{abOO\kO). 
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